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Abstract 

Background: Hunter syndrome (HS) is a lysosomal storage disease caused by iduronate-2-sulfatase (IDS) deficiency 
and loss of ability to break down and recycle the glycosaminoglycans, heparan and dermatan sulfate, leading to 
impairment of cellular processes and cell death. Cell activities and functioning of intracellular organelles are 
controlled by the clock genes (CGs), driving the rhythmic expression of clock controlled genes (CCGs). We aimed to 
evaluate the expression of CGs and downstream CCGs in HS, before and after enzyme replacement treatment with 
IDS. 

Methods: The expression levels of CGs and CCGs were evaluated by a whole transcriptome analysis through Next 
Generation Sequencing in normal primary human fibroblasts and fibroblasts of patients affected by HS before and 
24 h/144 h after IDS treatment. The time related expression of CGs after synchronization by serum shock was also 
evaluated by qRT-PCR before and after 24 hours of IDS treatment. 

Results: In HS fibroblasts we found altered expression of several CGs and CCGs, with dynamic changes 24 h and 
144 h after IDS treatment. A semantic hypergraph-based analysis highlighted five gene clusters significantly 
associated to important biological processes or pathways, and five genes, AHR, HIF1A, CRY1, ITGA5 and EIF2B3, 
proven to be central players in these pathways. After synchronization by serum shock and 24 h treatment with IDS 
the expression of ARNTL2 at 10 h (p = 0.036), PERI at 4 h (p = 0.019), PERI at 10 h (p = 0.041) and 16 h (p = 0.043) 
changed in HS fibroblasts. 

Conclusion: CG and CCG expression is altered in HS fibroblasts and IDS treatment determines dynamic 
modifications, suggesting a direct involvement of the CG machinery in the physiopathology of cellular 
derangements that characterize HS. 
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Background 

Hunter syndrome (HS), or Mucopolysaccharidosis type II, 
is a lysosomal storage disease and a X-linked recessive 
genetic disorder caused by a mutation in the IDS gene, 
leading to absence or deficiency of the enzyme iduronate- 
2-sulfatase (IDS). IDS deficit interferes with the ability to 
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break down and recycle glycosaminoglycans (GAGs), pre- 
cisely dermatan sulfate and heparan sulfate, causes lyso- 
somal engulfment, and hinders the customary working of 
cellular functions, compromising anatomical integrity in a 
number of organ systems and leading to curtailment of life 
expectancy [1]. The patients are characterized by coarsen- 
ing of facial features, bone and joint abnormalities, short 
stature, changes in the heart, respiratory system, hearing 
and vision, and in more severe forms by disturbed motor 
function, progressive learning difficulties and behav- 
ioural abnormalities [2]. Enzyme replacement therapy 
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with idursulfase, a recombinant form of human IDS, rep- 
resents the only affordable therapeutic approach presently. 

The cellular processes in every life form from bacteria 
to humans show diurnal variations driven by an internal 
timing system, and the oscillation frequency has a near- 
24-hour period, so that the rhythmicity is defined circa- 
dian (from the latin circa, about, and dies, a day) [3-5]. 
The mammalian circadian timing system is composed by 
a central pacemaker and master oscillator in the supra- 
chiasmatic nuclei (SCN) of the brain and self-sustained 
oscillators in the peripheral tissues [6,7]. The SCN are 
entrained to the environmental light/dark cycle by cues 
perceived through photon capture, signalling by 
melanopsin-containing retinal ganglion cells, and inputs 
conveyed by the retino-hypothalamic tract [8]. They 
synchronize peripheral oscillators and coordinate bodily 
functions by means of output pathways that may be 
neural projections (sympathetic and parasympathetic 
nerve fibers) or diffusible factors (for example hormones 
as Cortisol and melatonin) [9,10]. 

The molecular mechanisms underlying the biological 
clock functioning consist of a transcriptional-translational 
feedback loop operated by a set of genes, called core clock 
genes, coding for proteins that in turn block gene expres- 
sion, completing a cycle in approximately 24 hours [11]. 
The positive limb of the loop is represented by the 
CLOCK (or its paralog NPAS2) and ARNTL/BMAL1 (or 
its paralog ARNTL2/BMAL2) genes that code for the 
transcription factors CLOCK/NPAS2 and ARNTL/2. A 
heterodimer of the CLOCK/NPAS2 and ARNTL/2 pro- 
teins rhythmically activates the transcription of clock 
genes PERI, PER2, PER3, CRY1, CRY2, coding for PER 
and CRY proteins that represent the negative limb of 
the loop. PER and CRY proteins along with Casein 
Kinase IS/e (CKlS/e) and Casein Kinase II (CKII Al, 
A2, B), responsible of multilevel posttranslational regu- 
lation of various clock components, cooperate to form a 
repression complex that translocates back into the 
nucleus, interacts directly with CLOCK/NPAS2 and 
ARNTL/2, and impede their transcriptional activity 
[12]. The clock gene machinery is connected to a 
supporting feedback loop operated by the nuclear re- 
ceptors REV-ERBoc/[3 (encoded by the genes NR1D1 and 
NR1D2) and retinoic acid-related orphan receptor (ROR) 
a/y, where RORs positively regulate ARNTL/2 and REV- 
ERB expression and, in turn, REV-ERBs antagonize RORs 
[13,14]. Besides, SIRT1, a NAD + -dependent protein 
and histone deacetylase, is required for high-magnitude 
circadian transcription of several core clock genes, in- 
cluding ARNTL, PER2, and CRY1. SIRT1 counterbal- 
ances the histone acetyltransferase activity of CLOCK, 
binds CLOCK- ARNTL heterodimer in a circadian man- 
ner and promotes the deacetylation and degradation of 
PER2 [15]. 



TIMELESS is a core circadian clock gene in Drosophila 
melanogaster and is maintained in mammals, but its role 
in mammalian circadian clock function is not clear. 
TIMELESS and its partner TIMELESS interacting pro- 
tein (TIPIN) interact with components of the DNA rep- 
lication system to regulate DNA replication processes 
under both normal and stress conditions and are essential 
for ataxia telangiectasia and Rad3-related (ATR) -check- 
point kinase (Chk)l and ataxia telangiectasia mutated 
(ATM) -checkpoint kinase (Chk)2-mediated signaling and 
S-phase arrest [16-19]. 

The clock gene oscillation drives the rhythmic expression 
of other genes defined clock controlled genes, for example 
DBP, TEE, HLF, HSF1, NFIL3 (also called E4BP4), DEC1-2 
(also called BHLHE40-41), which steer downstream tissue 
specific genes regulating key cellular functions, such as cell 
cycle progression, proliferation, DNA damage response, au- 
tophagy, apoptosis, metabolism, redox equilibrium [20-23], 
and are involved in physiological processes such as inflam- 
mation/immune response [24,25], and paraphysiological 
phenomena, such as aging [26]. About 5% to 15% of 
genome-wide mRNA expression exhibits a circadian 
rhythm that is driven by the clock genes and there is a tis- 
sue specificity of cycling genes and a tissue specificity of 
clock controlled gene expression timing and level (output 
genes i.e., genes involved in the functional output of an 
organ) [27-30]. Accordingly, multicellular organisms are 
characterized by tissue-specific circadian regulation of tran- 
scription, with peripheral oscillators controlling the definite 
biochemical cascades relevant to their tissue or organ func- 
tion and generating rhythms in various pathways, including 
those involved in intracellular activities and depending on 
the correct functioning of intracellular organelles [27-30]. 

The importance of the role played by the clock gene ma- 
chinery in the regulation of the processes in the intracellu- 
lar organelles is confirmed by the physiopathology of 
Niemann-Pick types A and B disease, related to deficit of 
sphingomyelin synthase 2 (SGMS2), the enzyme that cata- 
lyzes the transfer of phosphocholine from phosphatidylcho- 
line onto ceramide to produce sphingomyelin, a major 
component of cell and Golgi membranes. Furthermore, 
Niemann-Pick type C disease is related to deficit of 
Niemann-Pick CI (NPC1), which is crucial for the intracel- 
lular trafficking of cholesterol from the late endosome to 
the trans-Golgi network [31]. SGMS2 and NPClprotein 
levels oscillate with circadian rhythmicity driven by the 
clock controlled genes SGMS2 and NPC1 respectively 
[27,28] . Patients carrying a mutation in these genes develop 
a condition characterized by accumulation of sphingomye- 
lin in spleen, liver, lungs, bone marrow and brain, causing 
irreversible neurological damage, and high cholesterol levels 
in the endosomal-lysosomal system, respectively [32] . 

The aim of our study was to assess the expression of clock 
genes and clock controlled genes in Mucopolysaccharidosis 
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type II, and to evaluate the circadian pattern of variation 
and the effects of the treatment with idursulfase on the ex- 
pression of clock genes and clock controlled genes at 
different time points. 

We addressed these issues taking advantage of an in vitro 
model of HS, represented by human fibroblasts with the 
mutational features of this mucopolysaccharidosis, through 
evaluation of the mRNA expression levels of the core clock 
genes and of a panel of clock controlled genes selected by 
means of literature-mining [27-30,33] (for a complete list of 
circadian transcripts refer to: CircaDB at http://circadb. 
org). These evaluations are part of a whole transcriptome 
analysis performed by Next Generation Sequencing (NGS) 
in fibroblasts from healthy subjects and fibroblasts from HS 
patients before, 24 hours and 144 hours after idursulfase 
treatment. We also evaluated by qRT-PCR the expression 
levels of core clock genes upon serum-shock induced 
synchronization in normal human fibroblasts and fibro- 
blasts of patients affected by HS before and after 24 hours 
of treatment with idursulfase. 

Materials and methods 

Cells 

Human fibroblasts from skin biopsy of five HS paediatric 
patients carrying different mutation in IDS gene were 
obtained from "Cell Line and DNA Bank from Patients 
Affected by Genetic Diseases", Gaslini Institute (Genoa, 
Italy). As healthy controls human fibroblasts from four 
children's circumcision were used; they were obtained 
from the Histology Unit of the Department of Histology, 
Microbiology and Medical Biotechnology (University of 
Padua, Padua, Italy). Written informed consents were 
obtained from patients at the time of biopsy and the 
study was approved by the Ethics Committee of the Uni- 
versity of Padua, Padua, Italy. All cells were anonym- 
ously obtained. 

Cell culture and treatment with idursulfase of HS 
fibroblasts 

Primary fibroblasts were cultured at 37° C in 5% C02 at- 
mosphere in RPMI medium supplemented with 15% 
fetal bovine serum (FBS), 100 U/ml penicillin and 
100 ng/ml streptomycin (Invitrogen Life Technologies, 
Milan, Italy). Fibroblasts of patients affected by HS were 
treated with idursulfase (Elaprase®, Shire Human Gen- 
etic Therapies, Inc, Lexington, MA, USA) at a concen- 
tration of 62.5 nM for 24 hours and the cells were 
harvested 24 h and 144 h after idursulfase treatment. 

Whole transcriptome analysis performed by NGS 
technology 

Total RNA was extracted using the TRIzol" Reagent 
(Sigma-Aldrich®) according to the manufacturer's proto- 
cols. Isolated RNA was quantified by NanoDrop ND- 



1000 Spectrophotometer (Thermo Scientific, Barnstead, 
NH, USA) and further assessed for quality using the 
Agilent 2100 Bioanalyzer (Santa Clara, CA) prior to li- 
brary construction. Total RNAs extracted from different 
cell lines were equally pooled into two specimens 
(Hunter and healthy) using 12.5 ug of RNA from each 
cell line. From total RNA of each pool mRNA purifica- 
tion was performed trough poliA(+) enrichment by 
Dynabeads® mRNA Purification Kit (Life Technologies'" 
Carlsbad, CA, USA). Total RNA was used for standard 
fragment-library preparation using the SOLiD Total 
RNA-Seq Kit (Life Technologies). Finally, emulsion PCR 
reactions were carried out by mixing appropriated 
amount of libraries with SOLiD beads. After PCR ampli- 
fication, emulsions were broken using butanol, and the 
beads were washed, enriched, and terminal transferased 
before quantification and deposition onto a slide for se- 
quencing. Templated beads were deposited onto one full 
slide, one sample for quarter (quad). Sequencing was 
carried out to 35 bases using SOLiD™ 3 (Sequencing by 
Oligo Ligation and Detection) System and following the 
manufacturer's instructions. 

Alignment 

Reads obtained from sequencing were aligned on human 
genome by CRIBI Biotechnology Centre (University of 
Padova) using the PASS software (http://pass.cribi.unipd.it). 
The February 2009 human genome assembly (GRCh37), 
publicly available at UCSC Genome Bioinformatics Site 
(http://genome.ucsc.edu/) was taken as reference sequence 
and the options best-hit, gap = 0 and maximun mismach = 3 
have been selected. Only unique reads (aligning on only 
one gene) and only genes with coverage higher than 50% 
of gene length were considered. 

Serum-shock induced synchronization and treatment with 
idursulfase of HS fibroblasts 

The serum shock induced synchronization was per- 
formed as follows: approximately 4x10 s cells/6 wells 
were plated the day before the experiments. Fibroblasts 
of patients affected by HS were treated with idursulfase 
at a concentration of 62.5 nM for 24 hours (Elaprase, 
Shire Human Genetic Therapies, Inc, Lexington, MA, 
USA). The day of the experiments, culture medium was 
exchanged with serum-rich medium (RPMI containing 
50% FBS) and after 2 hours this medium was replaced 
with serum free RPMI [34]. The cells were harvested 
1 h, 4 h, 10 h, 16 h, 22 h and 28 h after serum shock. 

Quantitative real time polymerase chain reaction 
(qRT-PCR) 

Total RNA was extracted from normal primary human fi- 
broblasts and fibroblasts of patients affected by HS before 
and after 24 hours of idursulfase treatment at the indicated 
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time points using the RNeasy" Mini Kit (Qiagen S.p.a. 
Milan, Italy) and subsequently digested by DNase I. Quanti- 
tative Real Time PCR was performed starting from 100 ng 
of purified RNA using the one step quantifast SYBR Green 
RT PCR KIT (Qiagen). For real-time PCR, we used the fol- 
lowing SYBR Green QuantiTect Primers purchased from 
Qiagen: ARNTL (QT00068250), ARNTL2 (QT00011844), 
CLOCK (QT00054481), CRY1 (QT00025067) PERI 
(QT00069265), PER2 (QT00011207), PER3 (QT00097713). 
Expression levels of target gene were normalized using the 
housekeeping control gene TATA binding protein (TBP, 
QT00000721). Values of mRNA expression levels of clock 
genes were calculated using the formula 2~ A . 

Statistical analysis 

Results are expressed as means ± SE of at least three differ- 
ent experiments. Comparisons were made using Students 
t-test as appropriate. The limit of statistical significance was 
set at to p < 0.001 for comparisons of levels of mRNA ex- 
pression determined by NGS, and p < 0.05 for comparisons 
of levels of mRNA expression determined by qRT-PCR. 
Each time series of clock gene expression levels was ana- 
lyzed for circadian rhythm characteristics by the single 
cosinor procedure involving the fit of a 24-h cosine curve 
to the data by least squares linear regression, in order to 
accurately describe waveforms and rhythm characteristics. 
An R2 value and a p value for the rejection of the zero- 
amplitude assumption were determined for each compo- 
nent in the cosine model separately and overall, with 
rhythm detection considered statistically significant if 
p < 0.05 and borderline significant if p > 0.05 and p < 0.10 
for any period tested. Circadian (24-h) characteristics were 
summarized from the 24-h cosine. Rhythm characteristics 
determined from the best-fitting cosine model include the 
MESOR (the middle of the cosine representing an adjusted 
average if unequal sampling), the amplitude (half the dis- 
tance from the peak and trough of the best fitting curve), 
and the acrophase of the cosine model (the peak of a single 
component cosine). All analyses were performed using 
the MATLAB 6.5 statistical package (MathWorks, Natick, 
MA, USA). 

Results 

The mRNA expression levels of the core clock genes 
ARNTL, ARNTL2, CLOCK, CRY1, CRY2, CSNK1D, 
CSNK1E, CSNK2A1, CSNK2A2, CSNK2B, NPAS2, NR1D1, 
NR1D2, PERI, PER2, PER3, RORA, SIRT1, TIMELESS, 
TIPIN and of a set of clock controlled genes were evaluated 
as part of a whole transcriptome analysis performed by 
NGS technology in healthy human fibroblasts (C) and 
fibroblasts of HS patients before (H), 24 hours (Tl) and 
144 hours (T2) after idursulfase treatment. Ratio, log2 
(ratio), fold change and p-values are shown in Additional 
file 1: Table SI. 



The expression levels of ARNTL, ARNTL2, CLOCK, 
CRY1, PERI, PER2, and PER3 were also evaluated by 
qRT-PCR after serum-shock induced synchronization in 
normal human fibroblasts and fibroblasts of patients af- 
fected by HS before and after 24 hours of treatment with 
idursulfase. 

Evaluation by NGS technology of core clock gene 
expression in HS fibroblasts vs healthy fibroblasts (H vs C) 

In HS fibroblasts the genes CLOCK (fold change = -1.45, 
p < 0.0001), NPAS2 (fold change = -1.65, p < 0.0001), 
CRY1 (fold change = -1.51, p < 0.0001), NR1D2 (fold 
change = -1.97, p < 0.0001) and SIRT1 (fold change = -1.30, 
p < 0.0001) showed lower expression levels compared to 
healthy controls. Conversely the genes PERI (fold change = 
1.27, p < 0.0001), CSNK1D (fold change = 1.26, p < 0.0001), 
CSNK1E (fold change = 1.54, p < 0.0001), and NR1D1 
(fold change = 1.99, p < 0.0001) showed higher expression 
levels in HS fibroblasts. The expression level of ARNTL 
(fold change = 1.34, p = 0.01), ARNTL2 (fold change = 
1.09, p = 0.47), PER2 (fold change = 1.40, p = 0.001), PER3 
(fold change = -1.01, p = 0.79), CRY2 (fold change = 1.16, 
p = 0.02), CSNK2A1 (fold change = 1.02, p = 0.56), 
CSNK2A2 (fold change = -1.08, p = 0.26), RORA (fold 
change = -1.07, p = 0.28), TIMELESS (fold change = 1.25, 
p = 0.006), and TIPIN (fold change = -1.25; p = 0.27) was 
not different in a statistically significant way between nor- 
mal and HS fibroblasts or did not reach the p < 0.001 
threshold value. 

Evaluation by NGS technology of core clock gene 
expression levels in HS fibroblasts after 24 hours of 
treatment with idursulfase versus healthy fibroblasts 
(Tl vs C) 

In HS fibroblasts after 24 hours of treatment with 
idursulfase NPAS2 (fold change = -1.39, p = 0.0005), 
CRY1 (fold change = -1.75, p < 0.0001) and NR1D2 (fold 
change = -1.65, p < 0.0001) genes showed lower expres- 
sion levels in comparison to controls. The genes ARNTL 
(fold change = 1.53, p < 0.0001), ARNTL2 (fold change = 
2.22, p < 0.0001), PER2 (fold change = 1.72, p < 0.0001), 
CSNK1E (fold change = 1.47, p < 0.0001), CSNK2B (fold 
change = 1.15, p < 0.0001), NR1D1 (fold change = 3.59, 
p < 0.0001), and TIMELESS (fold change = 2.41, p < 0.0001) 
showed higher expression levels in HS fibroblasts after 24 h 
treatment, in comparison to healthy human fibroblasts. 

The expression level of CLOCK (fold change = -1.13, 
p = 0.003), PERI (fold change = 1.22, p = 0.004), PER3 
(fold change = 1.14, p = 0.03), CRY2 (fold change = 1.04, 
p = 0.55), CSNK2A1 (fold change = 1.04, p = 0.30), 
CSNK1D (fold change = 1.12, p = 0.007), CSNK2A2 (fold 
change = -1.09, p = 0.25), RORA (fold change = -1.14, 
p = 0.03), SIRT1 (fold change = -1.12, p = 0.07) and TIPIN 
(fold change = 1.4; p = 0.03) was not different in a 
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statistically significant way between normal and 
idursulfase treated HS fibroblasts or did not reach the p < 
0.001 threshold value. 

Evaluation by NGS technology of core clock gene 
expression levels in HS fibroblasts 144 hours after 
treatment with idursulfase versus healthy fibroblasts 
(T2 vs C) 

In HS fibroblasts 144 hours after treatment with 
idursulfase the genes CLOCK (fold change = -1.66, 
p < 0.0001), NPAS2 (fold change = -1.56, p < 0.0001), 
CRY1 (fold change = -2.00, p < 0.0001), NR1D2 (fold 
change = -2.53, p < 0.0001), RORA (fold change = -1.28, 
p = 0.0003), TIMELESS (fold change = -1.93, 
p < 0.0001) and SIRT1 (fold change = -1.41, p < 0.0001) 
showed lower expression levels in comparison to controls. 
CSNK1E (fold change = 1.27, p < 0.0001), and NR1D1 
(fold change = 1.72, p < 0.0001) showed higher expression 
levels in HS fibroblasts 144 hours after idursulfase treat- 
ment in comparison to healthy human fibroblasts. The 
expression level of ARNTL (fold change = 1.38, p = 0.007), 
ARNTL2 (fold change = -1.36, p = 0.03), PERI (fold 
change = 1.02, p = 0.70), PER2 ((fold change = 1.20 
p = 0.10), PER3 (fold change = -1.08, p = 0.25), 
CRY2 (fold change = 1.05, p = 0.48), CSNK1D (fold 
change = 1.11, p = 0.01), CSNK2A1 (fold change = -1.14, 
p = 0.02), CSNK2A2 (fold change = -1.29, p = 0.002), 
CSNK2B (fold change = 1.00, p = 0.85), and TIPIN (fold 
change = -1.64; p = 0.03) was not different in a statisti- 
cally significant way between normal and idursulfase 
treated HS fibroblasts or did not reach the p < 0.001 
threshold value. 

Evaluation by NGS technology of changes of clock gene 
expression levels in HS fibroblasts after 24 hours of 
treatment with idursulfase versus untreated HS 
fibroblasts (T1 vs H) 

In HS fibroblasts after 24 hours of treatment with IDS 
CSNK2B (fold change = -1.60, p < 0.0001) showed lower 
expression levels compared to untreated HS fibroblasts 
(H). Conversely the genes CLOCK (fold change = 1.28, 
p < 0.0001), ARNTL2 (fold change = 2.03, p < 0.0001), 
NR1D1 (fold change = 1.80, p < 0.0001), NR1D2 (fold 
change = 1.19, p < 0.001), and TIMELESS (fold change = 
1.91, p < 0.0001) showed higher expression levels in HS 
fibroblasts treated for 24 hours with idursulfase , in 
comparison to untreated HS fibroblasts. The expression 
level of NPAS2 (fold change = 1.18, p = 0.11), ARNTL1 
(fold change = 1.14, p = 0.21), PERI (fold change = -1.04, 
p = 0.50), PER2 (fold change = 1.23, p = 0.02), PER3 (fold 
change = 1.16, p = 0.02), CRY1 (fold change = -1.16, p = 
0.09), CRY2 (fold change = -1.11, p = 0.11), CSNK1D (fold 
change = -1.12, p = 0.05), CSNK1E (fold change = -1.05, 
p = 0.20), CSNK2A1 (fold change = 1.01, p = 0.66), 



CSNK2A2 (fold change = -1.00, p = 0.98), RORA (fold 
change = -1.07, p = 0.30), SIRT1 (fold change = 1.16, 
p = 0.04), and TIPIN (fold change = 1.83; p = 0.01) 
was not different in a statistically significant way between 
treated and not treated HS fibroblasts or did not reach the 
p < 0.001 threshold value. 

Evaluation by NGS technology of changes of clock gene 
expression levels in HS fibroblasts after 144 hours of 
treatment with idursulfase versus untreated HS 
fibroblasts (T2 vs H) 

In HS fibroblasts after 144 hours of treatment with IDS 
the genes CSNK1E (fold change = -1.12, p < 0.0001), 
CSNK2A1 (fold change = -1.17, p < 0.001), CSNK2B 
(fold change = -1.21, p < 0.0001), NR1D2 (fold 
change = -1.28, p < 0.0001), and TIMELESS (fold 
change = -2.4, p < 0.0001) showed lower expression 
levels, in comparison to untreated HS fibroblasts. 

The expression level of CLOCK (fold change = -1.14, 
p = 0.009), NPAS2 (fold change = 1.05, p = 0.64), 
ARNTL (fold change = 1.03, p = 0.78), ARNTL2 (fold 
change = -1.48, p = 0.004), PERI (fold change = -1.24, 
p = 0.002), PER2 (fold change = -1.16, p = 0.14), PER3 
(fold change = -1.06, p = 0.37), CRY1 (fold change = -1.32, 
p = 0.002), CRY2 (fold change = -1.10, p = 0.15), 
CSNK1D (fold change = -1.13, p = 0.004), CSNK2A2 
(fold change = -1.18, p = 0.04), NR1D1 (fold change = -1.16, 
p = 0.01), RORA (fold change = -1.19, p = 0.01), SIRT1 (fold 
change = -1.08, p = 0.33), and TIPIN (fold change = -1.31; 
p = 0.26) was not different in a statistically significant way be- 
tween treated and not treated HS fibroblasts or did not reach 
the p < 0.001 threshold value. 



Evaluation by NGS technology of clock controlled gene 
expression in healthy fibroblasts, untreated HS fibroblasts 
and HS fibroblasts 24 hours and 144 hours after 
treatment with idursulfase 

Regarding to the clock controlled genes, in HS fibro- 
blasts there was a statistically significant difference and 
dynamic change with idursulfase treatment at the time 
points considered of the expression of a huge number of 
output genes involved in the control of important cellu- 
lar and tissue processes. The cellular processes driven by 
clock controlled genes differentially expressed with a 
statistical significance reaching the p < 0.001 threshold 
are represented by DNA transcription (AHR, ARNT, 
CTNNB1, FOS, FOXL2, FOXOl, HIF1A, HSF1, ID2, 
JUN, KEAP1, KITLG, KLF10, NFIL3, PPARA, PPARD, 
PPARG, RARA, RXRA, SMAD4, SP1, SREBF1, SRF, 
STAT1, STAT5A), post-translational modification and 
degradation (FBXL3, SUM03, USP5), biosynthesis 
(ACSL1, AEBP1, ALAS1, APBB1, APBB1IP, APBB2, 
BMPR1A, CEBPB, CREBBP, FASN, GYS1, HMGCR, 
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INSIG2, LPIN1, NAMPT), processing {ADA, ALDH1A1, 
ALDH1A3, ALDH1B1, ATF2, ATM, ATR, ATRIP, 
BHLHE40, BHLHE41, BNIP3, CES1, CES2, CYP7B1, E2F1, 
EGFR, EGR1, EIF2B3, HSP90AA1, HSPA1A, HSPAS, 
HSPD1, GABARAPL1, GL01, GLUT, HK1, LDHA, 
LMAN1, LMAN2, MAOA, NPC1, PARP1, PCK2, PRDX2, 
PYGL, SGMS2, ULK1, XBP1), transport {AP2A1, AP2M1, 
IGFBP3, IGFBPS, LDLR, SLC22A15, SLC25A1, SLC27A1, 
SLC7A8, SLC9A3R2, SLC9A9), DNA damage response 
(ATM, ATR, ATRIP), and cell cycle control (CCNA2, 
CCNB1, CCND1, CDK2AP1, CDKN1A, GADD4SA, 
MDM2,WEE1). The tissue processes influenced by 
deregulated clock controlled genes are mainly represented 
by inflammation, hemocoagulation and fibrinolysis 
(ADAM17, A2M, FN1, ICAM1, IL6, ITGA5, MASP1, 
MGST1, NFKBIA, PDGFRA, PDGFRB, PTGS2, SERPINE1, 
SPP1, THBD, TFPI2, TRAF2, VEGFA). 



Evaluation by qRT-PCR after Serum-Shock Induced 
Synchronization of changes of clock gene expression levels 
in healthy fibroblasts, untreated HS fibroblasts and HS 
fibroblasts after 24 hours of treatment with idursulfase 

After serum shock synchronization and qRT-PCR ana- 
lysis, statistically significant differences were evidenced 
in the mRNA expression levels of ARNTL2 at 10 h 
(p = 0.036), PERI at 4 h (p = 0.019), PER2 at 10 h 
(p = 0.041) and 16 h (p = 0.043), between idursulfase 
treated HS fibroblasts and untreated HS fibroblasts 
(Figure 1). Fitting of cosine curves with a 24-hour period 
to raw data and plotting as polarograms clock gene 
expression levels in control fibroblasts and fibroblasts of 
patients affected by HS before and after 24 hours of 
idursulfase treatment (Figures 2, 3 and 4) evidenced 
statistically significant and borderline significant rhythms 
of clock gene expression, and advance of phase os- 
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Figure 1 x-y plots representing the time related profiles of expression of ARNTL, ARNTL2, CLOCK, CRY1 PER1, PER2, and PER3 after 
synchronization with serum shock in control fibroblasts and fibroblasts of patients affected by Hunter Syndrome before and after 
treatment with iduronate-2-sulfatase. Original units standardized to TO (1 h after serum shock) and combined for analyses. * p < 0.05 for 
ANOVA followed by all pairwise multiple comparison procedures with Student-Newman-Keuls Method. 
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cillation after idursulfase treatment. MESOR, amplitude, 
acrophase and p values are reported in Table 1. 

Semantic hypergraph-based analysis of circadian gene 
expression 

The aim of this analysis was to determine the func- 
tional groups of genes that play any key-role in the 
development of HS. The analysis workflow is here 
preliminary sketched: (i) RNA-Sequencing: Detection 
of core clock gene and clock controlled gene expres- 
sion profiles; (ii) Descriptive analysis and feature se- 
lection: Selection of the most varying genes by the 
ensemble of principal component analysis and linear 
discriminant analysis; (iii) Qualitative hypergraph 
building: Inference of weighted and signed hypergraphs; 
(iv) Quantitative undirected graph construction: Math- 
ematical agglomeration of weights over the edges; (v) 
Weighted topology indices calculation: Computation 
of several weighted topology indices for each node 



and normalization of individual fold changes by the 
most discriminant index; (vi) Determination of the 
strongest connected components: Hierarchical cluster- 
ing and matching with critical biological functions 
[35,36]. 



Descriptive analysis and features selection 

We have preliminarily assessed the expression variation of 
clock genes and clock-controlled genes, through five differ- 
ent states: HS (H) vs. Control (C), Treatment at time 1 (Tl) 
vs. C, Treatment at time 2 (T2) vs. C, and Tl vs. H, T2 vs. 
H. We discarded genes that did not exhibit any relevant ex- 
pression change between these states. We combine a prin- 
cipal component analysis to determine those genes that 
maximally vary through the first three components, and a 
linear discriminant analysis to project the genes and their 
expressions onto the most varying components' axes 
(Figures 5 and 6). 



Control fibroblast 1 
Control fibroblast 2 




Figure 2 x-y plots showing from top to bottom the fitted cosine curves (continous line) superimposed on raw data (squares and 
circles) and polarograms of ARNTL, ARNTL2, CLOCK, CRY1 PERI, PER2, and PER3 expression levels after synchronization with serum 
shock in control fibroblasts. Original units standardized to TO (1 h after serum shock) and combined for analyses. Polarograms of cosinor 
analysis show the acrophases for the clock gene expression values. Radial axis represents the time point (in degrees) after serum shock 
corresponding to the acme and vector length represents the amplitude of the oscillation. 
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pjtient 1 ARNTL2 

Patient 1 + Elaprsse 




Figure 3 x-y plots showing from top to bottom the fitted cosine curves (continous line) superimposed on raw data (squares and 
circles) and polarograms of ARNTL, ARNTL2, CLOCK, CRY1 PER1, PER2, and PER3 expression levels after synchronization with serum 
shock in the fibroblasts of a patient affected by Hunter Syndrome before and after 24-h treatment with iduronate-2-sulfatase. Origina 
units standardized to TO (1 h after serum shock) and combined for analyses. Polarograms of cosinor analysis show the acrophases for the clock 
gene expression values. Radial axis represents the time point (in degrees) after serum shock corresponding to the acme and vector length 
represents the amplitude of the oscillation. 



Qualitative hypergraphs building 

We have assembled a hypergraph connecting genes by 
querying and merging a number of heterogeneous data 
sources: (i) Interpro and PFAM (protein domains), (ii) Gene 
Expression Omnibus (Co-localization and co-expression), 
(iii) BIOGRID and IREF (genetic interactions), (iv) 
PathwayCommons, IMID, NCI_NATURE, REACTOME, 
KEGG and BIOCARTA (pathways), (v) BIOGRID, BIND, 
HPRD, INTACT, MINT, MPPI and OPHID (physical inter- 
actions) and (VI) curated literature (predicted interactions). 
Two genes were connected by an edge, whenever at least 
an interaction evidence of any of the above mentioned 
interaction categories was found. Several pair of genes 
resulted to be connected by more than one edge. We built 
a weighted and signed hypergraph with weights over the 
edges (carrying the reliability of the corresponding interac- 
tions) and weights over the nodes (carrying the fold change 
expression of the corresponding genes) [37-39]. 



Quantitative undirected graph 

We deterministically inferred an undirectedm graph 
from our hypergraph, by applying an infective function 
to the sets of weights over the edges connecting any two 
pairs of genes. This simple agglomerative weighting 
function takes the weights of the edges connecting any 
two nodes in input, and gives a unique value in output. 
Constitutively, it gives more and more importance to the 
pairs that are connected by multiple edges. 

e" 

w AB = -y j w ABi 

I— I 

where n is the number of edges connecting any two 
nodes A and B, and i refers to the ith edge. 

Thus, we obtain a new graph that contains the same set 
of nodes, and individual edges connecting any two nodes, 
with the newly calculated weights over them (Figure 7). 
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Patient 2 ARNTL2 

Patient 2 ♦ Elaprase 




Figure 4 x-y plots showing from top to bottom the fitted cosine curves (continous line) superimposed on raw data (squares and 
circles) and polarograms of ARNTL, ARNTL2, CLOCK, CRY1 PER1, PER2, and PER3 expression levels after synchronization with serum 
shock in the fibroblasts of a patient affected by Hunter Syndrome before and after 24-h treatment with iduronate-2-sulfatase. Origina 
units standardized to TO (1 h after serum shock) and combined for analyses. Polarograms of cosinor analysis show the acrophases for the clock 
gene expression values. Radial axis represents the time point (in degrees) after serum shock corresponding to the acme and vector length 
represents the amplitude of the oscillation. 



Weighted topology indices calculation 

We calculated the closeness topology index for each 
gene. We recall that the un-weighted formulation of 
closeness between any two nodes of a graph measures 
the length of their shortest paths. The concept of close- 
ness originates from the inverse definition of farness, 
namely as the inverse of the sum of the distance from a 
node to all other nodes. Thus, the more important a 
node is, the lower is its total distance to all other nodes. 
Closeness can be regarded as a measure of how fast it 
will take to spread information from a node to all other 
nodes sequentially. Contrarily, the weighted formulation 
of closeness takes into account paths whose lengths are 
not the sum of the minimum number of hops from a 
node to all other nodes, but the sum of the inverse 
weights over the edges of the shortest paths from a node 
to all other nodes. This formulation gives importance to 
the reliability of a connection, when determining the 



essentiality of each node, thereby discriminating between 
reliable and not-reliable interaction paths. 

We have subjected the top 26 closest genes to other 
well-known topological indices (i.e., degree, betweenness, 
clustering coefficient and topological importance) and 
found that 5 (i.e., AHR, HIF1A, CRY1, ITGA5 and 
EIF2B3, see the plot below) were highly ranked by 4 out 
of 5 indices (Table 2 and Figure 8). 

Strongly connected components determination 

After having scaled the individual fold changes with the 
closeness values of each gene, we have finally determined 
the collaborative actions of such genes by extracting the 
most cohesive clusters out of their interaction network. We 
performed a hierarchical clustering on the expression pro- 
files of clock and clock controlled genes, by considering the 
Squared Euclidean distance as metric. We have subjected 
the found clusters to functional enrichment analysis against 
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Table 1 Parameters calculated by Cosinor analysis (fitting a cosine curve with a 24-hour period to raw data) of clock 
gene expression levels after serum-shock induced synchronization in normal fibroblasts and fibroblasts of patients 
affected by Hunter syndrome before and after 24 hours of treatment with idursulfase 



Fibroblast of control subject 1 






Fibroblast of control subject 2 








MESOR 


Amplitude 


Acrophase 


P 




MESOR 


Amplitude 


Acrophase 


P 


ARNTL 


0.578 


0.401 


60° 42' 


0.048 


ARNTL 


0.872 


0.095 


130° 30' 


0.747 


ARNTL2 


0.753 


0.186 


131° 58' 


0.734 


ARNTL2 


0.802 


0.041 


72° 33' 


0.948 


CLOCK 


0.834 


0.061 


45° 00' 


0.749 


CLOCK 


0.642 


0.096 


11° 11' 


0.827 


PERI 


0.813 


0.106 


51° 40' 


0.592 


PER1 


0.621 


0.108 


336° 40' 


0.818 


PER2 


0.777 


0.271 


277° 45' 


0.234 


PER2 


0.772 


0.030 


164° 21' 


0.983 


PER3 


0.420 


0.469 


258° 57' 


0.136 


PER3 


0.645 


0.185 


241° 35' 


0.825 


CRY1 


0.747 


0.184 


324° 45' 


0.577 


CRY1 


0.704 


0.091 


319° 46' 


0.842 


Fibroblast of Hunter syndrome patient 1 not treated 




Fibroblast of Hunter syndrome patient 1 treated with idursulfase 




MESOR 


Amplitude 


Acrophase 


P 




MESOR 


Amplitude 


Acrophase 


P 


ARNTL 


0.709 


0.190 


1 39° 04' 


0.493 


ARNTL 


0.499 


0.406 


69° 1 1 ' 


0.444 


ARNTL2 


0.466 


0.383 


204° 48' 


0.280 


ARNTL2 


0.597 


0.110 


101° 49' 


0.928 


CLOCK 


0.600 


0.218 


232° 42' 


0.825 


CLOCK 


0.500 


0.255 


37° 52' 


0.492 


PERI 


0.813 


0.207 


261° 52' 


0.516 


PER1 


0.321 


0.257 


28° 38' 


0.677 


PER2 


0.600 


0.282 


315° 00' 


0.294 


PER2 


0.585 


0.196 


330° 15' 


0.509 


PER3 


0.553 


0.193 


268° 32' 


0.867 


PER3 


0.484 


0.350 


305° 17' 


0.390 


CRY1 


0.461 


0.253 


161° 40' 


0.643 


CRY1 


0.399 


0.237 


2° 15' 


0.593 


Fibroblast of Hunter syndrome patient 2 not treated 




Fibroblast of Hunter syndrome patient 2 treated with idursulfase 




MESOR 


Amplitude 


Acrophase 


P 




MESOR 


Amplitude 


Acrophase 


P 


ARNTL 


0.552 


0.365 


103° 12' 


0.298 


ARNTL 


0.462 


0.407 


83° 41' 


0.163 


ARNTL2 


0.680 


0.310 


114° 09' 


0.260 


ARNTL2 


0.660 


0.219 


96° 41 ' 


0.349 


CLOCK 


0.662 


0.316 


119° 17' 


0.078 


CLOCK 


0.723 


0.235 


64° 40' 


0.296 


PERI 


0.597 


0.169 


54° 28' 


0.713 


PER1 


0.576 


0.323 


70° 31' 


0.235 


PER2 


0.420 


0.196 


9° 22' 


0.708 


PER2 


0.535 


0.229 


307° 31' 


0.682 


PER3 


0.383 


0.257 


322° 27' 


0.664 


PER3 


0.449 


0.308 


305° 32' 


0.572 


CRY1 


0.585 


0.502 


116° 45' 


0.009 


CRY1 


0.694 


0.302 


13° 30' 


0.039 



Cosinor, least-squares fit of single component cosine to all data. MESOR, rhythm-adjusted overall 24 h mean; Amplitude, peak-trough range of cosine model ; 
Acrophase, peak of cosine in degrees. 



the GO FAT category, and determined five clusters signifi- 
cantiy associated to at least a biological process or pathway 
(Figure 9). 

Cluster 1, which is made of CDK1, PER2, CCNA2, 
/UN, PERI, TIMELESS, CCNB1, CSNK2B, EIF2B3, 
NR1D1, ARNTL2, HSP90AA1, CSNK1D, E2F1, 
BHLHE40, CEBPB, TIPIN, ELAVL1, ATRIP, POLA1, 
MDM2, POLD1, HSPD1, HK1, LDHA, TRAF2, HSPA5, 
TBP, PPP5C, ALDH1A3, CCRN4L, CDC7, SUM03, 
AHCY, POR, GNB2L1, ESRRA, ATP6V1D, USP5, 
DDX11, NRF1, FOXL2, is strongly associated with sev- 
eral cellular processes: cell cycle (GO:0022402, p-value = 
1.46E-02; GO:0000278, p-value = 2.40E-03; etc.), DNA 
damage response activities (GO:0000077, p-value = 
1.055E-3; GO:0007050, p-value = 4.634E-3, etc.) and 
pathways, such as cell cycle (WikiPathways, WP179, 
p-value = 2.443E-4), and hypoxia-inducible factor in 



the cardiovascular system (BIOCARTA, p-value = 
2.041E-2). 

Cluster 2, which is made of ARNTL, RXRA, CREBBP, 
ITGA5, SERPINE1, CCND1, CSNK1E, BHLHE41, 
CDKN1A, FN1, PKM2, PTGS2, DBP, PPARD, HSF1, 
THBD, IGFBP3, HSPA1A, AP2M1, TUBB, FASN, 
P4HA2, ALAS1, ALDH1B1, CDK2AP1, GYS1, NOTCH1, 
ADA, SLC7A8, SCG5, is associated to the inflammatory 
response (GO:0006954, p-value = 3.043E-2), response 
to lipid (GO:0033993, p-value = 1.346E-3), liver devel- 
opment (GO:0001889, 1.025E-3) and cell motility 
(GO:0048870, p-value = 3.889E-3; GO:0016477, p-value = 
1.823E-3). 

Cluster 3, which is made of WEE1, FOS, EGR1, 
PDGFRB, FOXOl, PPARG, IL6, NFYB, A2M, PDGFRA, 
PTGDS, ICAM1, IGFBP5, NFIL3, AEBP1, LGALS9, is in- 
volved in the response to cold (GO:0009409, p-value = 
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Figure 5 Heat map of the core clock gene expression fold changes between untreated patients and control cases (H/C), patients 24 or 
144 hours after treatment and control cases (T1/C and T2/C), treated and untreated patients (T1/H and T2/H). Cyan segments and a 
density plot quantify the magnitude of intra class and global FCs, respectively. 



3.126E-5) and glucose homeostasis (GO:0042593, NFYA, RORA, CHEK2, ATM, NONO, KLF10, RBX1, 

p-value = 5.443E-3). TPMT, KITLG, TTC1, SEC61G, SEC61A1, GCLC, 

Cluster 4, which is made of CLOCK, ARNT, PER3, SLC29A1, RCC1, SLC9A8, APBB2, SRP9, CHMP1B, 

EIF2BS, NPAS2, CSNK2A1, CSNK2A2, GSK3B, ATR, ENOX2, SLC9A9, ITFG2, WDR5, AACS, cooperate with 

HDAC3, SIRT1, EIF2B2, PLAT, PARP1, HSF2, MYC, genes in cluster 1 to negatively regulate cell cycle 



Mazzoccoli et al. BMC Medical Genomics 201 3, 6:37 
http://www.biomedcentral.eom/1755-8794/6/37 



Page 12 of 20 






Li 



CD? r.i ; 

TTGA5 

p.- 

: E 
F7G.II 

= :>>r ■ 



NOTCH 1 
El_i -T. " 
HK1 



. iSPOi 

I C 

I f:-_a- 

I cc--- 




I ALDHTA3 
I 

I BMLHE+0 



■ , - : - 

IP ="3 

PRKCA 
E'.tF"' - 
Hfinl 
RB1 
GLUL 



ACSL1 
OCL-2 
NPC1 
GL"'* 

ME 1 1 

= - - f :, 

HAHFT 

SLCft*S 
A2M 
Cur 
E'.p- 



H/C Tl/C T2/C Tl/H T2/H 



Figure 6 Heat map of the clock-controlled gene expression fold 
changes between untreated patients and control cases (H/C), 
patients 24 or 144 hours after treatment and control cases (T1/ 
C and T2/C), treated and untreated patients (T1/H and T2/H). 

Cyan segments and a density plot quantify the magnitude of intra 
class and global FCs, respectively. 



(GO:0007050, p-value = 3.094E-3), to the response to 
decreased oxygen levels (GO:0036293, p-value = 
3.128E-2) and to participate to the P53 signaling 
pathway (Pathway Ontology, PW:0000718, p-value = 
1.466E-2). 

Cluster 5, which is made of AHR, HIF1A, CRY1, RBI, 
STAT1, NFKBIA, STAT5A, PRKCA, FBXL3, NR1D2, 
NFE2L2, HES1, ACSL1, CXCL12, NNMT, NAMPT, 
EGFR, FMR1, NPC1, PYGL, NACA, LDLR, GLUL, 
GLOl, MGST1, BMPR1A, plays an important role in the 
positive regulation of epithelial cell proliferation 
(GO:0050678, p-value = 8.024E-7), response to lipid 
(GO:0033993, p-value = 3.898E-4; GO:0071396, p-value = 
3.709E-3), tissue morphogenesis (GO:0048729, p-value = 
7487E-4) and takes part in the adipogenesis pathway 
(WikiPathways, WP236, p-value = 8.700E-4). 

Discussion 

The molecular oscillator ticking in every cell regulates 
the customary sequence of intracellular activities, 
allowing the coordination of key pathways and the 
compartmentalization in the temporal dimension of 
poorly-compatible biochemical processes [40]. Altered 
functioning of the clock gene machinery may determine 
loss of time-of-day specific transcription of clock genes 
and clock controlled genes, causing severe deregulation 
of cellular homeostasis, cell dysfunction and biochemical 
and structural derangements that may lead to cell death 
and tissue dysfunction [41]. 

The aim of our study was to evaluate the functioning 
of the clock gene machinery in HS, and to address this 
issue we analyzed and compared mRNA expression 
levels of the core clock genes CLOCK, NPAS2, ARNTL1, 
ARNTL2, PERI, PER2, PER3, CRY1, CRY2, CSNK1D, 
CSNK1E, CSNK2A1, CSNK2A2.CSNK2B, NR1D1, 
NR1D2, RORA, SIRT1, TIMELESS, TIPIN, and a set of 
clock controlled genes in normal human fibroblasts 
and fibroblasts derived from subjects affected by 
Mucopolysaccharidosis type II. We also evaluated 
in vitro the effects of the treatment with idursulfase on 
circadian gene expression at different time points. 

Analyzing data obtained from NGS we observed al- 
tered mRNA levels of some of the clock genes and clock 
controlled genes and a variable response to idursulfase 
treatment. 

In HS fibroblasts CLOCK and NPAS2 showed lower 
expression levels. CLOCK and its paralog NPAS2 encode 
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(See figure on previous page.) 

Figure 7 Undirected and weighted interaction network between core clock genes and clock controlled genes in five different states: 
Hunter disease (H) vs. Control (C) (H/C), Treatment at time 1 (T1) vs. C (T1/C), Treatment at time 2 (T2) vs. C (T2/C) and T1 vs. H (T1/H), 
T2 vs. H (T2/H). Nodes colors range from red (negative FC) to green (positive FC). Diameter of nodes is proportional to the magnitude of the FC. 
Thickness of edges varies with the reliability of the interaction. Overall, the network exhibits the following global topological properties: clustering 
coefficient = 0.293, connected components = 1, diameter = 4, radius = 2, centralization = 0.259, shortest paths = 35910 (100%), characteristic 
path length = 2.002, average number of neighbors = 25.484, density = 0.135 and heterogeneity = 0.498. 



transcription factors belonging to the basic helix-loop- 
helix/Per-Arnt-Simpleminded (bHLH-PAS) domain fam- 
ily, and represent key elements in the positive limb of 
the transcriptional-translational feed-back loop that 
drives circadian rhythmicity of the oscillatory functions 
underlying cell homeostasis. In particular, their encoded 
proteins have histone acetyltransferase activity and are 
involved in chromatin remodeling with specificity for 
histones H3/H4 [42]. ARNTL, the heterodimerization 
partner of CLOCK and NPAS2, enhances their enzym- 
atic function, and the decreased expression of these cen- 
tral components of the molecular clock could hinder the 
correct clock gene expression in the fibroblasts of 
Mucopolysaccharidosis type II patients. 

A lower expression level was found in HS fibroblasts for 
CRY1, whereas PERI was up-regulated. The mammalian 
Period and Cryptochrome genes are both E-box-regulated 
genes, but in peripheral tissues CRY1 mRNA expres- 
sion peak is delayed by several hours with respect to 
that of PERI. CRY1 usually shows evening-time ex- 
pression and serving as a strong repressor at morning- 
time elements ensures a delay in feedback repression 
in the molecular clock. This phase delay in CRY1 tran- 
scription is required for customary mammalian clock 
function [43]. The E-box (CACGTG) is crucial for day- 
time transcriptional activity and the delay originates 
from interactions between the proximal E-box and 
ROR response elements (RORE) present in the CRY1 
promoter [44]. 

A higher expression level was found in HS fibroblasts 
for CSNK1D and CSNK1E, encoding the serine/threo- 
nine protein kinase casein kinase (CK) I 5 and e respect- 
ively, which are key regulators of metazoan circadian 
rhythmicity. CKI binds to and phosphorylates the Period 
proteins, wich in turn interact with a variety of circadian 
regulators, suggesting the possibility that CKI may inter- 
act with and phosphorylate additional clock components 
as well [45]. Cryptochrome proteins are phosphorylated 
by CKI only when both proteins are bound to mamma- 
lian Period proteins [12]. ARNTL is also a substrate for 
CKI in vitro, and CKI activity positively regulates 
ARNTL-dependent transcription from circadian pro- 
moters in reporter assays [12]. CKI phosphorylates mul- 
tiple circadian substrates and may exert its effects on 
circadian rhythm in part by a direct effect on ARNTL- 
dependent transcription [12,45]. The up-regulation of 



these protein kinases found in HS fibroblasts could play 
a crucial role in the derangement of the clock gene ma- 
chinery and the downstream controlled pathways, but 
on the other hand CKs could represent important 
molecular targets for new therapeutic approaches. 

In HS fibroblasts we found higher expression of 
NR1D1, and lower expression of NR1D2. NR1D1 and 
NR1D2, activated by CLOCK/NPAS2ARNTL/ARNTL2 
heterodimers, code respectively for the orfan nuclear re- 
ceptors REV ERBa and REV ERBp 1 , which realize a supple- 
mentary loop amplifying and stabilizing the oscillation of 
the molecular clockwork. REV-ERBa and REV-ERB|3 act 
as negative transcriptional regulators by binding RORE 
in gene promoters, preventing the binding of the posi- 
tive transcriptional regulator RORct, and negatively 
regulating the expression of ARNTL, CLOCK, and CRY1 
[46,47]. In addition, REV-ERBa is directly involved in 
lipid metabolism antagonizing the opposite role of 
RORa and inhibiting the expression of genes coding for 
apolipoprotein C-III, a constituent of very-low-density 
lipoproteins [48]. 

A lower expression level was found in HS fibroblasts 
for SLRT1, which regulates metabolic and stress re- 
sponses acting in concert with the circadian rhythm ma- 
chinery [49]. SIRT1 is involved in a number of cellular 
processes, including gene silencing at telomere and mat- 
ing loci, and modulates cell survival by inhibiting apop- 
tosis or cellular senescence induced by external 
challenges, including DNA damage and oxidative stress 
[50]. SIRT1 impinges on Ku70, peroxisome proliferator 
activated receptors, p53 and the forkhead box O 
(FOXO) family of transcription factors, and is modulated 
by active regulator of SIRT1 (AROS), hypermethylated 
in cancer 1 (HIC-1), deleted in breast cancer 1 (DBC1) 
and E2F transcription factor 1 (E2F1) [49-51]. Interest- 
ingly, E2F1 was greatly over-expressed in HS fibroblasts, 
but normalized after 144 hours of IDS treatment, 
whereas FOXOl was severely down regulated, and 
responded faintly to IDS treatment. Besides, in HS fibro- 
blasts there was severe down-regulation with no re- 
sponse to IDS treatment of NAMPT, encoding the 
enzyme controlling the synthesis of nicotinamide aden- 
ine dinucleotide (NAD + ), a cofactor of SIRT1, linking 
nutrient sensing and circadian regulation [52]. 

After 24 hours of treatment with idursulfase in fibro- 
blasts of patients with Mucopolysaccharidosis type II the 



Table 2 The top 26 closest genes subjected to other well-known topological indices (i.e., degree, betweenness, clustering coefficient and topological 
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Figure 8 Genes with the highest topological scores. AHR, HIF1A, CRY1, EIF2B3, and ITGA5 are among the top 26 most ranked genes by 4 out 
of 5 indices. 



expression evaluated by NGS of the genes CLOCK, 
ARNTL2, NR1D1, NR1D2, and TIMELESS showed 
higher expression levels compared to untreated HS fi- 
broblasts. This effect faded completely after 144 hours, 
and at this time point CSNK1E, CSNK2A1, CSNK2B, 
NR1D2, and TIMELESS genes showed decreasing ex- 
pression levels. Accordingly, the evaluation by qRT-PCR 
of clock gene expression upon synchronization through 
serum shock of normal fibroblasts and fibroblasts of 



patients affected by Hunter syndrome before and after 
24 hours of idursulfase treatment evidenced a statisti- 
cally significant effect on the expression levels of 
ARNTL2, PERI and PER2, but only at the earlier time 
points examined, with advance of the phase of oscilla- 
tion of gene expression. 

Regarding to the clock controlled genes, in Hunter syn- 
drome fibroblasts there was alteration and dynamic change 
with idursulfase treatment at the time point considered of 




Figure 9 Hierarchical clustering of the weighted interaction networks. Genes are clustered according to the proximity of their expression 
profiles, evaluated by an Euclidean distance metrics. 
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the expression of a huge number of output genes involved 
in the control of important cellular and tissue processes. 
An example is represented by SERPINE1, which showed 
higher expression levels in HS fibroblasts when compared 
to healthy fibroblasts, increased after 24 hours of 
idursulfase treatment, and decreased after 144 hours. 
SERPINE1 codes for the plasminogen activator inhibitor 
(PAI)-l, the major component of inhibitors of fibrinolysis, 
whose activity shows a clear circadian oscillation peaking in 
the early morning. CLOCK:ARNTL and CLOCK:ARNTL2 
heterodimers activate the SERPINE1 promoter, driving the 
circadian variation in circulating PAI-1, and CLOCK: 
ARNTL2 heterodimer, binding two E-box enhancers in 
the promoter, shows double capability to activate PAI-1 
expression [53]. 

The timing of processes that underlay cell function and 
the customary sequence of activation of key pathways de- 
pend on the correct ticking of the biological clock. Accord- 
ingly, the altered expression of clock controlled genes 
found in HS fibroblasts could hinder normal output of the 
clock gene machinery {BHLHE40, BHLHE41, NFIL3), and 
cellular activities, such as DNA transcription (CTNNB1, 
FOS, FOXL2, FOXOl, HIF1A, HSF1, ID2, JUN, KEAP1, 
KITLG, KLF10, PPARA, PPARD, PPARG, RARA, RXRA, 
SMAD4, SP1, SREBF1, SRF, STAT1, STATSA), post- 
translational modification and degradation (FBXL3, 
SUM03, USP5), lipid and glucose metabolism and biosyn- 
thetic pathways (ACSL1, AEBP1, ALAS1, APBB1, APBB1IP, 
APBB2, BMPR1A, CREBBP, FASN, GYS1, HMGCR, 
INSIG2, LPIN1, NAMPT), molecular processing {ADA, 
ALDH1A1, ALDH1A3, ALDH1B1, ATF2, CES1, CES2, 
CYP7B1, E2F1, EGFR, EGR1, EIF2B3, HSP90AA1, 
HSPA1A, HSPA5, HSPD1, GLOl, GLUL, HK1, LDHA, 
LMAN1, LMAN2, MAOA, NPC1, PARP1, PCK2, PRDX2, 
PYGL, SGMS2), molecular transport [AP2A1, AP2M1, 
IGFBP3, IGFBP5, LDLR, SLC22A15, SLC25A1, SLC27A1, 
SLC7A8, SLC9A3R2, SLC9A9), DNA damage response 
(ATM, ATR, ATRIP), endoplasmic reticulum stress and un- 
folded protein response (XBP1), xenobiotic response (AHR, 
ARNT), autophagy (BNIP3, CEBPB, GABARAPL1, ULK1), 
cell cycle control (CCNA2, CCNB1, CCND1, CDK2AP1, 
CDKN1A, GADD45A, MDM2, WEE1), and tissue pro- 
cesses, such as inflammation, hemocoagulation and fibrin- 
olysis (ADAM17, A2M, FN1, ICAM1, IL6, ITGA5, MASP1, 
MGST1, NFKBIA, PDGFRA, PDGFRB, PTGS2, SERPINE1, 
SPP1, TFPI2, THBD, TRAF2,VEGFA) [27-30,54]. 

The semantic hypergraph-based analysis of circadian 
gene expression brought out five gene clusters significantly 
associated to at least one biological process or pathway (cell 
cycle, DNA damage response, inflammation, liver develop- 
ment, cell motility, glucose homeostasis, response to cold 
and to decreased oxygen levels, the P53 signaling pathway, 
the positive regulation of epithelial cell proliferation, tissue 
morphogenesis, response to lipid, and the adipogenesis 



pathway). Besides, five genes, AHR, HIF1A, CRY1, ITGA5, 
and EIF2B3, were highlighted as top ranked by all the five 
considered topological indices and characterized by multi- 
faceted centrality nature in the network analysis of the cir- 
cadian transcriptome. The results suggest an interesting 
connection between circadian and hypoxic pathways and 
the toxicological signal mediator, the Aryl Hydrocarbon Re- 
ceptor (AHR). Upon xenobiotic binding, AHR transcrip- 
tionally activates xenobiotic metabolizing enzymes and 
regulates molecules involved in the signaling of nuclear 
factor-erythroid 2-related factor-2 (Nrf2), p53 (TRP53), ret- 
inoblastoma (RBI), and NF-kB, influencing cellular re- 
sponses to oxidative stress and inflammation. Through 
activation of these signaling pathways, regulation of cell 
cycle, and interactions with hypoxia-inducible factors 
(HIFs), AHR takes part in endogenous developmental func- 
tions during hematopoietic stem-cell maintenance and 
differentiation [55]. Intriguingly, the deficit of CRY pro- 
teins causes constitutive NF-kB and protein kinase A sig- 
naling activation, leading to unrelenting increase of 
proinflammatory cytokines, and ITGAS encodes a receptor 
for fibronectin and fibrinogen, playing a role in chemotaxis, 
leukocyte migration, cell-cell adhesion, angiogenesis, and 
blood coagulation [56,57]. 

On the other hand, a cross-talk exists between circadian 
and AHR signaling pathways at genetic, epigenetic and 
proteomic level, and this interaction might play a role in 
the regulatory influences maneuvered by circadian rhythmi- 
city on cell physiology [58]. Regarding to the interaction be- 
tween hypoxic pathways and AHR complex, it has been 
proven that in endothelial cells AHR has a physiological 
function in the absence of exogenous ligands, and upon ac- 
tivation mediates toxicity of halogenated aromatic hydro- 
carbons. In addition, evidence suggests that hypoxia 
induces HIF dependent gene expression, significantly redu- 
cing AHR expression level and AHR-mediated expression 
of cytochrome P-450 phase I xenobiotic metabolizing en- 
zyme CYP1A1 [59]. A cross-talk exists between hypoxic 
and circadian pathways operated by the PAS protein family 
members PER and CLOCK and HIF-la. Interestingly, the 
components of AHR, circadian and hypoxic pathways are 
characterized by a PAS domain that serves as an interface 
for protein-protein interactions [60]. Moreover, cell prolif- 
eration and differentiation both in vitro and in vivo are 
highly influenced by oxygen concentration, and neuronal 
stem cell proliferation and neuronal and oligodendroglial 
differentiation are enhanced by a mild level of hypoxia [61]. 
Intriguingly, EIF2B3 encodes one of the subunits of initi- 
ation factor EIF2B, which catalyzes the exchange of 
eukaryotic initiation factor 2-bound guanosine-diphosphate 
(GDP) for guanosine-5'-triphosphate (GTP), representing 
an essential factor for protein synthesis, and mutations in 
this gene have been associated with neurodegenerative and 
white matter diseases [62-64]. 
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Talcing into account some limitations of the study, repre- 
sented by the limited number of patients, suffering from a 
disease that is clinically heterogeneous in terms of onset, se- 
verity and progression, the dynamic variation of expression 
levels observed for core clock genes and clock controlled 
genes in the fibroblasts of patients affected by HS at differ- 
ent time points of treatment with idursulfase may be an in- 
direct evidence of the key role played by the molecular 
clock in the regulation of the complex array of cellular 
functions in this mucopolysaccharidosis. A paradigm is rep- 
resented by the altered expression of the clock controlled 
genes NPC1 and SGMS2, which are down-regulated in 
HS fibroblasts and increase temporarily after idursulfase 
treatment. 

Conclusion 

Biological processes show time related variations driven 
by genetically encoded oscillators operated by transcrip- 
tional/translational feedback loops hardwired by the 
clock genes and their coded circadian proteins. Data 
obtained in our study show that circadian gene expres- 
sion is variably altered in the fibroblasts of subjects af- 
fected by HS, the treatment with idursulfase determines 
some modifications in the expression levels of clock 
genes and clock controlled genes, but these changes are 
temporary and fade in the course of treatment. Altered 
functioning of the clock gene machinery in HS may de- 
termine severe deregulation of cellular homeostasis, cell 
dysfunction and biochemical and structural derange- 
ments that may lead to cell death. These alterations are 
related to the key role played by the molecular clockwork 
in the control of downstream gene expression regulating a 
complex array of cellular functions, such as molecule bio- 
synthesis, post-translational modification, processing, 
transport, conjugation, internalization and degradation, 
and cell processes such as cell cycle, autophagy, apoptosis 
and DNA damage response.The alteration of clock gene 
expression levels and the response to the deficient enzyme 
suggest a direct involvement of the molecular clock ma- 
chinery in the physiopathology of cellular derangements 
that characterize Mucopolysaccharidosis type II, opening 
the way to possible new therapeutic strategies. 
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Additional file 1 : Table SI . Expression levels of core clock genes and 
clock controlled genes evaluated by whole transcriptome analysis 
performed through Next Generation Sequencing technology in normal 
human fibroblasts (Control) and fibroblasts of mucopolysaccharidosis 
Type II patients before (Hunter), 24 hours (T1) and 144 hours (T2) after 
idursulfase treatment. 
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